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For a second-order phase transition the critical energy range of interest is larger than the energy 
range covered by a canonical Monte Carlo simulation at the critical temperature. Such an extended 
energy range can be covered by performing a Wang-Landau recursion for the spectral density followed 
by a multicanonical simulation with fixed weights. But in the conventional approach one loses the 
advantage due to cluster algorithms. A cluster version of the Wang-Landau recursion together 
with a subsequent multibondic simulation improves for 2D and 3D Ising models the efficiency of 
the conventional Wang-Landau/multicanonical approach by power laws in the lattice size. In our 
simulations real gains in CPU time reach two orders of magnitude. 
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Equilibrium properties of statistical physics systems 
are often estimated by Markov chain Monte Carlo 
(MCMC) simulations Q. In many cases one is interested 
in calculating expectation values for a range of temper- 
atures with respect to the Gibbs canonical ensemble. It 
has turned out that instead of performing simulations 
of the canonical ensemble at distinct temperatures it is 
often advantageous to combine them into one simula- 
tion of a "generalized" ensemble @, H, 0, [f| ; for reviews 
see dHll. 

While the power of generalized ensembles is well docu- 
mented for first-order phase transitions and complex sys- 
tems such as spin glasses and peptides (small proteins), 
this is not the case for second-order phase transitions, 
although convenience of such applications is claimed by 
Landau and collaborators [9( . However they lose the cru- 
cial advantage which cluster algorithms [Ifl, El provide 
for MCMC simulations of second-order phase transitions. 
Here we present a generalization to cluster algorithms. 
To keep the paper accessible for non-experts, we restrict 
our investigations to 2D and 3D Ising models, while the 
points made are generally valid for cluster algorithms. 

In MCMC simulations of second-order phase transi- 
tions one wants to cover the scaling region in which many 
physical observables diverge with increasing lattice size. 
So we ask the question: How large is the energy range of 
this region on a finite lattice and is it eventually already 
covered by a single canonical simulations at the (infinite 
volume) critical temperature T c = 1//3 C ? 

For simplicity our lattices are of shape L D and periodic 
boundary conditions are assumed. We denote the prob- 
ability density of the energy from a canonical MCMC 
simulation by P(E). Finite-size scaling (FSS) arguments 
[l2| imply C ~ L a l v for the specific heat at /3 C , where 
a and v are, respectively, the critical exponents of the 
specific heat C and the correlation length £. A second- 
order phase transition requires v > 0. Let us first assume 



a > 0. The fluctuation-dissipation theorem gives 

(E - E) 2 ^ ~ L d+*/v where E = (E), (1) 

implying for the range covered by the simulation at (3 C 

AE = \E Q . 75 - E . 25 1 ~ L D l 2+a ' 2v , (2) 

where E q , q = 0.25 and q = 0.75, are fractiles of the 
energy distribution [Tj. In the vicinity of (3 C (A constant) 



E(P)/L D = E{f3 c )/L D +A(j3 — (3 c f 



(3) 



and using the hyperscaling relation [12j a — 2 — Dv, 
Eq. ||3J| translates into a reweighting range 



A/3 ~ L- 1 '" 



(4) 



The desired reweighting range is determined by the need 
to cover the maxima of all divergent observables mea- 
sured. Let the maximum value of such an observable 
S L (0) be S^ ax = Sl(Pl ax ) and denote the critical expo- 
nent of S by a. Then FSS theory implies 



gmax 



L a/u 



(5) 



Reweighting has to cover a reasonable range about the 
maximum, say from f3 L ~ to /3£ + > f3 L ~ defined as solu- 
tions of 



S L ((3) = rSr x , 0<r<l. 



(6) 



which becomes large for r small. We define (3 r L € 
{/3£~,/3£ + } to be the (3 r L value which is further away 
from [3 C than the other and assume 



A/?l = 1/31 -P c \=a r L-\ 



(7) 



where a r and k > are constants (k independent of ?*). 
For sufficiently large L we suppose that 



SlWD = S reg + A (A/3£)" 



(8) 
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FIG. 1: Canonical (indicated by "rwght") versus desired (en- 
tire f3 axis) reweighting range on an 80 3 lattice. 



holds, where S rcs is a regular background term. Com- 
bining Eqs. |5]), ([7|), and ((8|) we conclude 



(9) 



i.e., the desired range ([7]) scales in the same way as the 
canonical range (Tjf]) . However, the proportionality factor 
a r can be much larger than the one encountered for the 
canonical range. With the modest value r = 2/3 this 
point is made in Fig. [JJ for the 3D Ising model on an 80 3 
lattice. We plot the specific heat C(0) and for S((3) the 
first structure factor (see, e.g., Ref. [l3|), whose maxi- 
mum scales ~ L 1 I V . The desired reweighting range is 
more than 17 times larger than the canonical reweight- 
ing range from a simulation at /3™ ax of the specific heat 
(in realistic applications one does not know /3 C a-priori 
and /?™ ax of a suitable observable is a good substitute). 

Using the same line of arguments for a logarithmic sin- 
gularity 



S((3) = S^-A\n\/3-(3 c \ 



(10) 



one finds that the exponent n in Eq. ([7]) is no longer 
independent of r, but 



-/u 



(11) 



While the canonical reweighting range scales still ~ 
L- 1 ^, the desired reweighting range becomes ~ L r ' v ', 
so that the ratio desired/canonical diverges ~ jj. 1 - r )l v . 
With S = C the 2D Ising model provides an example. 

In conclusion many more than one canonical simu- 
lation are typically needed to cover a relevant part of 
the scaling region of a second-order phase transition. In 
principle this can be done by patching canonical sim- 
ulations from several temperatures together, relying on 
a multi- histogram approach [lij ]. Besides that dealing 
with many simulations is tedious, weaknesses of these 
approaches are that the histograms fluctuate indepen- 
dently and that their patching has to be done in regions 



where the statistics is reduced due to the decline of the 
number of histograms entries. More stable estimates are 
obtained by constructing a generalized ensemble, which 
allows the random walker to cover the entire region of 
interest. This requires two steps: 

1. Obtain a working estimate of the weight factors. 

2. Perform a MCMC simulation with fixed weights. 

To be definite we confine our discussion to the multi- 
canonical (MUCA) simulations Extension to cluster 
algorithms are known (l5l . [111. We will rely on multi- 
bondic (MUBO) simulations T5[. This defines step 2, 
but leaves still many options open to deal with step 1. 
"Working estimate" means that the approximation of the 
weights of the generalized ensemble is good enough so 
that the energy range in question is covered in step 2. 
Historically step 1 has been one of the stumbling blocks 
of umbrella sampling: "The difficulty of finding such 
weighting factors has prevented wide applications of the 
umbrella sampling method to many physical systems" 
Most convenient is to have an efficient general pur- 
pose recursion for step 1 at hand. While designs were re- 
ported in a number of papers see also Refs. 0,11, Hfl, 
the winning approach appears to be the one of Wang 
and Landau (WL) Q (although somewhat surprisingly 
we found only one comparative study [l9j). 

The WL recursion differs fundamentally from the ear- 
lier approaches by iterating the weight at energy E multi- 
plicatively with a factor /wl > 1 rather than additively. 
At a first glance the WL approach is counter-intuitive, 
because the correct iteration of the weight factor close to 
the desired fixed point is obviously proportional to one 
over the number of histogram entries H(E) and not to 

l//wi , ■ However, what matters is a rapid approach 
to a working estimate. The advantage of the WL over 
the other recursions is that it moves right away rapidly 
through the targeted energy range. When it is sufficiently 
covered, the iteration factor is refined by /wl ^ V/wl, 
so that it approaches 1 rapidly. Once the system cycles 
with frozen weights through the desired energy range the 
goal of a working estimate has been reached and the WL 
recursion is no longer needed (20| . We now generalize 
this approach to cluster algorithms. 

We use the energy function of the g-state Potts models, 



£ = -2£<W 

(ij) 



(12) 



where the sum is over the nearest-neighbor sites of a D- 
dimensional cubic lattice of N = L D Potts spins, which 
take the values qi — 1, . . . ,q. The factor of two has been 
introduced so that q = 2 matches with the energy and [3 
conventions of the Ising model literature. 

In the Fortuin-Kasteleyn (FK) cluster language (2l| 
the Potts model partition is written as 



]T z({ii}AK}) with 



3 



n [ as "iu 



+ <5, 



& y oJ 



(13) 



where a = e 2/3 — 1. For a fixed configuration {qi} of Potts 
states the Swendsen-Wang updating procedure is to 
generate bonds variables bij (simply called bonds in the 
following) on links with 8 qiqj — 1: Bonds bij = 
generated with probability p and bonds bij ■ 
probability q so that p/q — a and p + q = 1 holds, 
gives p 



1 - e" 2/3 for fe, : 
0. On <L 



1 and q = 1 — p 
links we have b' S4 



1 are 
with 
This 

e -2/3 

for 6^ = 0. On = links we have b'^ = with 

probability one. We call bonds with bij = 1 active or set. 
A cluster of spins is defined as a set of spins connected 
by active bonds and an update is to flip entire clusters of 
spins, {qi} -> {<£}. 

Let us denote the number of active bonds by P. The 
MUBO partition function [l5[ is defined by 

Zmubo -EE W{ - B ) ( 14 ) 

{«}{&«} 

where a bond weight factor W(B) has been introduced. 
A valid updating procedure for the configurations of this 
partition function is formulated in the following. 

A. For qi ^ qj a bond is never set. This applies to 
the initial as well as to the updated bond on this link, 
so that B does not change. B. For qi = qj there are two 
possibilities: 



TABLE I: 3D Ising model simulations on L lattices. 



01 



a WL recursion production 



20 0.210 649 0.233 690 2" 18 
30 0.216 443 0.229 336 2~ 18 
44 0.218 545 0.227 013 2 -19 
56 0.219 755 0.225 914 2" 19 
66 0.220 063 0.224 709 2 -21 
80 0.220 482 0.224 377 2 -21 



19 962 
27 344 
33 266 
56 323 
62 884 
108 618 



32 x 32 768 
32 x 32 768 
32 x 65 536 
32 x 65 536 
32 x 131 072 
36 x 131 072 
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FIG. 2: ri n t(L) for the 3D Ising model (see text). 



1. The initial bond is not set, fey = 0. Then B' = B 



for ti 4 = and B' = B 



1 for fe'„. = 1. The 



updating probabilities are 



Pi(0^0) = 



q W(B) 



qW{B) +pW(B + l) 
and Pi (0^ 1) = 1-Pi(0 ->{)). 



(15) 



1. Then B' = B — 1 
for b'ij = and B' = B for b' tJ = 1. The updating 



2. The initial bond is set, fey 



probabilities are 
ft(l-»0) = 



qW(B - 1) 



qW(B - l)+pW(B) 
and P 2 (l -»■ 1) = 1 - P 2 (l -► 0). 



(16) 



After the configuration is partitioned into clusters [22| , 
the update is completed by assigning with uniform prob- 
ability a spin in the range 1, . . . , q to each cluster. 

In its generalization to cluster algorithms the WL re- 
cursion updates then In W(B) according to 

InW(B) -» lnW(P)-a W L, a WL = M/wl) , (17) 

whenever a configuration with B bonds is visited. All 
recursions are started with o^vl = 1 and we iterate 
a WL — * owl/2 according to the following criteria: 



1. The Markov chain just cycled from B\ to and 
back. Here ~W L and B r L are bond estimates corre- 
sponding to (3 r ^~ and respectively, determined 
by short canonical simulations. 

2. The bond histogram h(B), measured since the last 
iteration, fulfilled a flatness criterion /i m in//imax > 
cut, where cut was equal to 1 /3 in most of our runs. 

3. We freeze the weights after a last iteration is per- 
formed with the desired minimum value a^J. 



After a short equilibration run, measurements are 
performed during the subsequent simulation with fixed 
weights, each tuned to approximately 1 000 cycling 
events. Canonical expectation values at inverse tempera- 
ture (3, < (3 < /3£ + are obtained by reweighting (frj|. 
Table[T]gives an overview of our 3D Ising model statistics. 
The effectiveness of the recursion is seen from the fact 
that it takes never more than 3% percent of the statis- 
tics used for production (these numbers are in sweeps). 

Similarly the initial simulations, which determine B L 
and P^ , take less than 3%. 

From the production statistics we calculate integrated 
autocorrelation times r; n t and compare them in Fig. [2] 
with those of a MUCA simulation. From the MUBO time 
series we calculated r int for (a) energies and (b) bonds 
and found the results almost identical (slightly higher 
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FIG. 3: Tmt(-L) for the 2D Ising model. 



for the energies, but indistinguishable on the scale of the 
figure). For MUCA the estimates are from energies. Up 
to a constant factor practically identical results are ob- 
tained from cycling times. In our code one MUCA sweep 
was about three times faster than one MUBO sweep. 

The critical slowing down is ~ L z . For the dynamical 
critical exponent we find z = 2.22 (11) for MUCA and 
z = 1.05 (5) for MUBO. So the performance gain is a bit 
better than linear with the lattice size L. The data in 
Fig. [2] scatter more than one might have expected about 
the fits because our fi r £~ and /3£ + values are based on 
MCMC estimates, which are by themselves noisy. Our 
exponent for cluster updating is significantly higher than 
the one estimated from simulations at /3 C , z — 0.50 (3), for 



the Swendsen-Wang algorithm [23j. The reason is that 
the efficiency of the cluster algorithm deteriorates off the 
critical point, even when one is still in the scaling region. 
Therefore, we think that our exponent of z k, 1 reflects 
the slowing down in real application more accurately than 
the small value of the literature. In particular the cluster 
algorithm becomes rather inefficient for calculating the 
long tail of the specific heat for j3 > /3™ ax . 

In Fig. [3] we show integrated autocorrelation times 
from simulations of the 2D Ising model for which we ad- 
justed our simulation parameter to cover the full width 
at half- maximum of the specific heat. This corresponds 
to r = 1/2 in Eq. (fTTj) . The dynamical critical expo- 
nent takes than the values z — 2.50 (4) for MUCA and 
z = 1.04(2) for MUBO. The MUCA value reflects that 
the number of canonical simulations needed to cover the 
desired energy range grows now ~ L 1 / 2 , while the canon- 
ical critical value is approximately two @, [HJ . 

Finally we remark that the efficiency of simulations of 
second-order phase transitions can presumably be fur- 
ther improved by optimizing the weights with respect to 
cycling along the lines introduced in Ref. [25j. 
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